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Abstract. A 2D contact dynamics model is proposed as a microscopic description of a collapsing suspen- 
sion/soil to capture the essential physical processes underlying the dynamics of generation and collapse 
of the system. Our physical model is compared with real data obtained from in situ measurements 
performed with a natural collapsing/suspension soil. We show that the shear strength behavior of our 
collapsing suspension/soil model is very similar to the behavior of this collapsing suspension soil, for both 
the unperturbed and the perturbed phases of the material. 
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1 Introduction 
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In nature, fragile or metastable granular structures, col- 
lapsing under an applied load exist in a large variety. For 
instance, these collapsing suspensions/soils are well known 
from soil mechanics as an unsaturated (by water) grain 
structure which collapses when adding water (thus becom- 
ing saturated) or/and when loaded. They are also known 
to be loose structures and can be water deposited or ae- 
olian 1,2,3J. Similar structures are also responsible for 
collapsible behavior of colloidal gels [4] or even snow [5]. 

In the present work we develop and numerically solve 
a simple physical model for a collapsing suspension/soil 
representation, that is capable of incorporating at the mi- 
croscopic level the essential structural and dynamical fea- 
tures of the material when subjected to different types of 
perturbations. Within this context we focus on two as- 
pects of the problem, namely the penetration of an in- 
truder, leading to a partial collapse of the system, and 
the shear resistance property of the material. The model 
results of our approach are then compared to real and in 
situ data of a collapsing suspension/soil. 



2 Simulation method 

The discrete element method constitutes a general class of 
modeling techniques to simulate the microscopic behavior 
(i.e., at the particle scale) of granular/soil materials. Here 
a variant of the contact dynamics method, originally de- 
veloped to model compact and dry systems with lasting 
contacts [6], is used to describe a 2D cohesive suspen- 
sion/soil material. This technique is conceptually based 
on the exact implementation of non-smooth contact laws, 
which means that the steric volume exclusion for perfectly 



rigid particles and the Coulomb friction law are strictly 
implemented [IllH]. 

The absence of cohesion between particles can only be 
justified in dry systems on scales where the cohesive force 
is weak compared to the gravitational force on the parti- 
cle, i.e. for dry sand and coarser materials, which can lead 
to densities close to that of random dense packings. For 
collapsing suspensions/soils, an attractive force must play 
an important role in the stabilization of large voids [9j,fT0j , 
leading to highly porous systems as e.g. in fine cohesive 
powders. In the nanometer range of particle sizes, the co- 
hesive force becomes the dominant force, so that particles 
stick together upon first contact. Here the bonding be- 
tween two particles is considered in terms of a cohesion 
model with a constant attractive force F c acting within 
a finite range d c , so that for the opening of a contact a 
finite energy barrier F c d c must be overcome. The utiliza- 
tion of this simplified cohesion model is justified in the 
case of short range interactions as we are not interested in 
the detailed behavior of an individual contact, but on the 
collective behavior of the system, where the opening of a 
contact is mostly determined by the force and energy bar- 
rier. More precisely, using a more complex model would 
not lead to a drastic change in the macroscopic behavior 
when keeping the same force and energy barrier. In most 
force-driven simulation as in this paper the cohesion force 
dominates the behavior and the energy barrier plays a mi- 
nor role but its existence is of crucial importance to avoid 
simulation artefacts [SUHJ- As we show later, this cohesion 
force can be mediated, for example, by bacteria present 
in the system which are much smaller than the grains, 
i.e. we have an effective short range interaction. We use 
d c = 10 -4 (in units of particle radii) for all the follow- 
ing simulations. In addition, we implement friction and 
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rolling friction between two particles in contact, so that 
large pores can be stable in the system. For comparison, 
we adopt the same friction coefficient /x = 0.3 and coef- 
ficient for rolling friction \i r = 0.1 (in units of particle 
radii) in all simulations, taking into account that rolling 
of particle contacts is easier than sliding [9l[8] . 

In the case of collapsing suspensions/soils one also has 
to take into account the time necessary for bonds to ap- 
pear, i.e. during relatively fast processes new bonds will 
not be formed, whereas for longterm processes bonds are 
allowed to form at a particle contact. Finally, gravity also 
cannot be neglected in the model since the particle di- 
ameter is usually well above the micron-size. For simplic- 
ity, however, the surrounding pore water is not explic- 
itly considered but only taken into account as a buoy- 
ant medium, reducing the effective gravity acting onto 
the grains. In this way, we disregard fluid motion. This 
hypothesis should remain valid for the settling grains in 
the system as long as the velocities reached by each parti- 
cle are sufficiently large (i.e., inertial effects are relevant), 
and the total virtual buoyancy induced by all particles re- 
leased during the penetration process is small [H]. The 
issue of disregarding the fluid motion will be discussed in 
more detail when presenting the results in the following 
section. 



3 Results 

As mentioned earlier the major feature of the metastable 
granular structures, investigated in this paper is that they 
collapse, e.g. when an external load is applied. This will be 
illustrated and investigated in section 13.11 by penetration 
of an intruder leading to partial collapse of the structures. 
Understanding the penetration behavior will also be im- 
portant for the experiments presented in section [3?3l where 
the measurement device has to be pushed inside the soil 
before the measurements. Before and after the collapse 
the shear strength of the material changes drastically as 
shown in the simulations (sec. 13. 2\ and in experiments 

(sec, org} . 

A useful number for characterization of the structures 
after the (partial) collapse where many bonds break is 
the fraction of still active bonds (i.e. fully cohesive) with 
respect to the total number of contacts. This number will 
be used in the following for interpretation of the results. 



3.1 Penetration of an intruder 

We carry out simulations of a large disk (intruder) of low 
density (half the grain density) pushed with constant force 
into the granular structure and subsequently removed again 
from a fragile cohesive granular structure. Only if the force 
exceeds a certain threshold it penetrates the medium. Ad- 
ditionally, we assume that the time scale for the pack- 
ing generation is large enough to create bonds between 
the grains, whereas the movement of the intruder is on 




Fig. 1. Snapshots from simulations showing a typical real- 
ization of the penetration process of an intruder into a very 
loose cohesive packing. The (unperturbed) collapsing suspen- 
sion/soil shown in (a) is modeled as a tenuous granular net- 
work of cohesive disks assembled by a gravity driven process 
of ballistic deposition and contact dynamics [6]|9] . The system 
size is 51 x 199 in units of grain diameters, corresponding to 
1.7cmx6.7cm when using the typical diameter d = 344/xm of 
the grains in the experiments. As shown in (b), the movement 
of the intruder is responsible for the partial destruction of the 
granular structure along its trajectory. At the end (c), the in- 
truder rests under a compact mass of (perturbed) material. 



time scales much smaller than the bonding time. The se- 
quence of snapshots shown in Fig. [T] corresponds to differ- 
ent stages of the penetration process into a typical real- 
ization of a very loose packing. The initial (unperturbed) 
network structure (Fig.QJi) is a result of the aggregation of 
cohesive disks by means of a gravity driven process which 
includes ballistic deposition with contact dynamics [6,9J. 
At this point we are mostly interested in the penetration 
process into a fragile structure independent on its process 
of generation. The specific configurations used here are 
generated to resemble as closely as possible the experi- 
ments described later (sec. I3.3[) . The details of this gener- 
ation process will be described later (sec. 13. 2\i . The initial 
configurations used here (Fig. [1]) are created with a value 
Fc/ F g = 10 4 for the cohesion force F c normalized by the 
gravitational force F g on each grain. Deliberately we chose 
this high value to obtain a rather loose structure with 
coordination number slightly above 2 (2.06), since lower 
values of F c /F g lead to denser structures (cf. sec. l3.2p with 
smaller pores, i.e. a more homogeneous microstructure as 
found for non-cohesive materials and therefore relatively 
dense materials in laboratory experiments j 1 211 1 311 1 1| . 

Note, that the microstructure of two dimensional sys- 
tems gives, at first glance, a less homogeneous impres- 
sion than the corresponding three dimensional systems as 
pores are much easier to identify visually. This has been 
illustrated, e.g. for the compaction of cohesive powders 
when comparing two dimensional simulations [S] to three 
dimensional simulations [15] . 
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Fig. 2. Dependence of the normalized weight of material 
above the intruder after penetration w/wt (open circles, left 
ordinate) on the applied pushing force Fd/F c (normalized by 
cohesive force). In the same plot, we also show how the nor- 
malized force necessary to pull out the intruder F u /F c (filled 
circles, right ordinate) varies with Fd/F c . 



we show the dependence of the weight fraction w/wt on 
the ratio Fd/F c , where w is the weight of material above 
the intruder after the penetration process, wt is the to- 
tal weight of the sample, Fa is the force applied to push 
the disk down and F c is the characteristic value used for 
the cohesion force. As depicted, the force to pull up the 
intruder disk can indeed be as much as three times higher 
than the force to push it down. Note, that the fraction of 
active bonds in the final packing strongly depends on the 
"penetration strength" . For low values of Fd/F c , where the 
intruder is not pushed into the structure very deeply, most 
of the bonds are still active. For Fd/F c ~ 0.33, e.g. 93% 
of the bonds are active with respect to the final struc- 
ture, corresponding to 97% of the original bonds of the 
initial structure, i.e. very few bonds are broken. Also, for 
Fd/F c ~ 1, where the intruder almost reaches the bot- 
tom, still 68% of the bonds are active (78% of the original 
bonds). As the intruder can only go down until reaching 
the bottom, increasing the "penetration strength" further 
only leads to the breaking of more cohesive bonds. This in- 
dicates a stronger collapse of the structures which explains 
why the the force to pull up the intruder still increases for 
high values of Fd/F c . 



Once the intruder breaks through, the highly porous 
material collapses under the action of gravity. We observe 
the creation of a channel (Fig. \Vp) which finally collapses 
over the descending intruder. At the end, the partial de- 
struction of the network due to the penetration process 
generates a large amount of compact (perturbed) mate- 
rial over the intruder. As shown in Fig.QJ, the larger disk 
is buried under the debris of particles. In this case of a 
relatively "strong penetration" (cf. following paragraph) 
the fraction of active bonds for the total packing is about 
30%. This means that about 43% of the cohesive contacts 
in the initial configuration remain intact as during the 
collapse new contacts are formed, leading to an increase 
in the average coordination number to about 2.82. These 
numbers, averaged over the whole systems, can be seen as 
a measure for the extent of still intact areas and collapsed 
areas. In the fully intact areas still all bonds are active and 
the coordination number is unchanged, slightly above 2. 
On the contrary, in the fully collapsed area, i.e. the chan- 
nel created by the intruder, only about 2% of the bonds 
are still active, and the coordination number is about 3.2. 
In this respect the total numbers can also be viewed as an 
indicator for the collapsed fraction of the structure. 

The results presented above clearly suggests that ob- 
jects lighter than water can be effectively swallowed in 
a collapsing suspension/soil. Furthermore, if we assume 
that cohesive bonds can be restored by some particu- 
lar physico-chemical or biological mechanism, the force 
needed to remove the intruder disk could be significantly 
higher than the originally penetration force. In order to 
provide a more quantitative confirmation for this interest- 
ing behavior, several simulations of the intruder penetra- 
tion/removal process have been performed for 10 different 
realizations of granular networks generated with the same 
microscopic cohesion and friction parameters. In Fig. [2] 



The results presented above illustrate the typical be- 
havior of such a fragile structure for a given initial density 
(here the initial volume fraction is about 0.43). Chang- 
ing the volume fraction (density), the curves are shifted 
to the right for higher densities or to the left for lower 
densities, still showing the same qualitative behavior. We 
expect three dimensional simulations to display a quali- 
tative similar penetration behavior, but the intruder ex- 
periencing less resistance of the material due to the addi- 
tional degree of freedom. This effect has been found when 
comparing the pore stabilization mechanisms of two di- 
mensional and three dimensional simulations [9"lll5j. Dis- 
regarding the interstitial fluid motion keeps the model as 
simple as possible and nevertheless able to reproduce the 
main experimental observations. Of course the details of 
the collapse of the material may be influenced by the flow 
field of the surrounding fluid (e.g. [T4|fT3]). We tested the 
influence of the fluid by introducing a viscous drag on the 
grains. For this, the drag coefficient of water and a typical 
grain size of d = 344 /xm have been used. The simulation 
results showed no significant difference. 

The results of this section are also very important 
for the experiments presented in section T3.3I using a vane 
rheometer pushed into the investigated soil. When pushed 
in smoothly the thin vane blades will not damage too much 
the structure for the material to be considered unper- 
turbed. This is one reason for choosing a relatively small 
intruder, being only three times larger than the grains, 
and being in the range or even smaller than the typical 
pore size of the structure. Additionally, this saves a lot of 
computation time as for larger intruders also the system 
sizes have to be larger. Preliminary studies with larger in- 
truders show similar qualitative behavior and will be the 
subject of future studies. 
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3.2 Shear Resistance 

The unperturbed material is modeled by a ballistic depo- 
sition of particles driven by gravity, as an example of a 
process that generates very fragile structures. During this 
relatively slow process, bonds are allowed to form when 
particles stick to each other. As previously discussed the 
origin of the attractive force can be explained in terms 
of cohesive bonds mediated by the bacteria living in the 
suspension [16) . For the bonds, we use the same model as 
described before, i.e. the bond strength is F c , which is the 
dominant parameter determining the density of the sys- 
tem. The particles are sequentially deposited with a fixed 
time interval between each one to allow for relaxation. 
To get a very loose and fragile structure, an extremely 
large time between two depositions should be used. The 
procedure adopted here is to gradually increase the time 
interval. We only stop when the packing generated with 
the highest time interval has, within a given tolerance, 
the same particle density of the previous one. We show in 
Fig. [3^ that this procedure results in a highly porous, and 
therefore tenuously connected, network of grains. In prac- 
tice, for the simulations presented here, depositing one 
particle each time step turned out to be already adiabatic, 
i.e. the density did not change any more (less than 0.1%). 
Note, that the generation process of the fragile structure 
used in our simulations somehow resembles the processes 
of formation of natural collapsing suspensions / soils inves- 
tigated in in-situ experiments [16] . This will later be im- 
portant to be able to compare the experiments with our 
model (sec. I3.3[) . 

With the purpose to simulate the shear strength af- 
ter the collapse, we implement in our model a piston-like 
pushing force acting on the particles at the top. In this 
situation, the collapse of the granular network that starts 
from the upper region of the system is analogous to a 
gradual increase in the effective piston weight, leading to 
an acceleration of the compaction mechanism. During this 
process bonds will break without reforming again up to the 
point in which most of the bonds are broken (see Fig.[3jD). 
This behavior will be valid for all materials where the 
time scale of bond formation is larger than that of the 
collapsing and measurement, or where the bonds are not 
able to reform again at all, as e.g. in the in situ exper- 
iments (sec. 13 .3[) . After the collapse only 8.5% of the 
bonds are active (i.e. cohesive contacts) corresponding to 
about 13% of the active bonds of the initial structure, i.e. 
most of the bonds are broken. As the number of contacts 
increases during the collapse the coordination number in- 
creases from slightly above 2 (2.05) to about 3.17, whereas 
the volume fraction increases from about 0.4 to about 0.77. 
The value for the coordination number is almost the same 
as within the channel after the compaction by an intruder 
(cf. sec. 13.11 coordination number about 3.2), whereas the 
fraction of active bonds, although very low, is still higher 
than in the intruder case. This can be understood because 
the intruder destroys the bonds very efficiently on its way 
down. Then, the material settles again without almost no 
active bonds in the "channel" . In the case of compacting 
the whole system "externally" as described above not all 




Fig. 3. Typical realizations of the granular system pre- 
pared for computational simulations of shear strength (periodic 
boundaries in horizontal direction). In (a) we show the (unper- 
turbed) granular network of cohesive disks (F c /F g — 10 4 ) gen- 
erated by ballistic deposition and contact dynamics, and in (b) 
the (perturbed) compact structure obtained after the collapse 
of the system due to a piston-like pushing force applied at the 
top The system size is 51 x 199 in units of grain diameters, cor- 
responding to 1. 7cm x 6. 7cm when using the typical diameter 
d = 344/^m of the grains in the experiments, b: 51 x 106. 



bonds have to be broken to compact the system. For the 
case of very large intruders one expects both situtions: 
One part of the system below the intruder is compacted 
similar as the whole system. Another part of the system 
collapses over the intruder within its "channel" . 

We calculate the shear strength of the unperturbed 
and perturbed structure as follows. At a given depth, we 
apply a constant force in the horizontal direction to a ran- 
domly chosen particle and observe how far it can move in 
this direction. Our assumption here is based on the fact 
that in real experiments the rheometer creates a thin shear 
layer at the upper and lower edges of the vane. This layer 
is represented in our model by these sample particles sub- 
jected to a constant force. By changing the force using 
nested intervals, one can calculate within a given numer- 
ical tolerance the minimum force necessary to move this 
particle at a distance that is sufficiently large to sample 
the disordered porous geometry (e.g., approximately 20 
particle diameters away) . This procedure is then repeated 
for different particles at the same depth to produce an av- 
erage shear strength value. During this procedure, mea- 
suring the threshold in a static system when motion sets 
in, it is justified to disregard the fluid motion. 
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Fig. 4. Dependence of the simulated shear strength F s (nor- 
malized by the gravitional force F g of each grain) on the depth 
for unperturbed and perturbed systems (depth h in units of 
grain diameters d). Compared to the perturbed system (dia- 
monds), the unperturbed grain structure shows only a weaker 
dependence on the depth. This is similar to the behavior of 
a real collapsing suspension/soil (see Fig. [5]). For the unper- 
turbed system, an increase by a factor of 100 in the cohesion 
force does not change substantially the behavior of the shear 
strength. For the perturbed system we also show two different 
values of the friction coefficient jj, (0.3 and 0.6). A higher fric- 
tion coefficient leads to a stronger increasing shear strength as 
indicated by the slopes (a « 23.4 for fj, — 0.3 and a » 35 for 
H = 0.6) of the linear fit as in Eq.lfljl, 



The results in Fig. |4] show the variation with depth 
of the shear strength averaged over 10 different particles 
for both unperturbed and perturbed systems. The cohe- 
sion forces are dominant in the unperturbed state while 
the contribution of friction forces governs the dynamical 
behavior of the perturbed material. The model results 
obtained with the unperturbed system show only a minor 
dependence of the shear strength on the depth compared 
to the perturbed structure. Furthermore, for the unper- 
turbed system the shear strength does not extrapolate to 
zero in the limit of low depths. Additionally, we observe 
that a significant increase in the cohesion force intensity 
(e.g., by a factor of one hundred) does not lead to a sub- 
stantially different behavior of the shear strength. This is 
a consequence of two distinct effects, namely that (?) the 
particle density decreases with the cohesion force, and that 
(ii) a larger shear strength is necessary to move a particle 
in a more compact medium (i.e., a high density packing 
with volume fraction of about 0.69 and coordination num- 
ber about 3). As also shown in Fig. [4] the shear strength 
increases linearly with depth for the case of the perturbed 
material. Expectedly, as we observe in the same plot, an 
increase in the friction coefficient generally increases the 
average shear strength. Using the slopes of Fig. [4] one can 
for the perturbed case calculate the effective friction coeffi- 
cient of the material by simply deviding the shear strength 



at a certain depth by the weight of the particles above that 
depth. Here we use the full depth of the whole structure to 
achieve highest accuracy. For our simulation this results in 
an effective friction coefficient of 0.32 (fi = 0.3) and 0.48 
(/i = 0.6). The relatively high values and in particular the 
large difference can not solely be explained by the pres- 
ence of friction. It can better be understood by the fact 
that there are still active bonds present in the structure, 
enhancing the strength of the material. 

3.3 Comparison to Experiments 

In order to test our physical model we use experimental 
data from investigations of a specific type of a natural 
collapsing suspension/soil [IS]. This comparison between 
the real data and our contact dynamics model can only be 
made qualitatively due to obvious difficulties in obtaining 
reasonably precise experimental measurement of any in 
situ microscopic parameters. 

The experimental data come from a natural reserve 
called Lengois Maranhenses located in the North-East of 
Brazil [15]. We found that at the shore of many dry- 
ing lakes in this place, it is common to find a special 
type of a collapsing suspension/soil consisting of an im- 
permeable crust lying above a metastable suspension of 
grains [17]. The slow drying process of the lakes, includ- 
ing continuous deposition of grains transported by the 
wind from the adjacant dunes into the lakes, is consistent 
with the slow limit of depositing grains used in the simula- 
tions (sec. l3.2p leading to very fragile structures. Provided 
one does not exert on the surface a pressure higher than 
p c = 10 — 20 kPa, it is possible to step on it and the surface 
will elastically deform in a very similar way to what hap- 
pens when one walks on a waterbed. These deformations 
visibly extend over a couple of meters. If the pressure p c 
is exceeded, the surface cracks in a brittle way produc- 
ing a network of tensile cracks [16] . After that we observe 
the separation of the excess pore-water from a repacked 
and wet sand soil 18J with pronounced shear thinning be- 
havior [19]. Because the collapse of this suspension/soil 
is irreversible, we had to study its rheology and strength 
in situ. This irreversible collapse leading to a compacted 
material with reduced pore volume, could be qualitatively 
reproduced in our simulations, as presented before (cf. 
sec. O and EH]). 

By placing light plates on the surface we could walk 
on the suspension/soil without visually modifying it and 
make various measurements before and after the collapse. 
The most striking result concerns the shear strength r 
measured with a vane rheometer. The vane consists of 4 
thin blades (see inset in fig. [5]) . During the measurement 
the vane rotates cutting out a cylindrical area of the ma- 
terial. In this respect the results of section 13.11 are very 
important, as the rheometer has to be pushed into the 
investigated soil, before measuring. The thin vane blades 
ensure that when pushed in smoothly enough most of the 
structure in the cylinder rotated is not not destroyed by 
the vane. Thus, we were really able to measure the prop- 
erties of the unperturbed material. 
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Fig. 5. Experimental measurements of the shear strength 
as a function of depth before (open circles) and after (filled 
circles) the collapse of the suspension/soil. For the measure- 
ments, we used a vane rheometer Geonor H-60 shown in the 
inset. The horizontal error bars indicate the read-off accuracy 
of the vane rheometer, the horizontal error bars are given by 
the height of the vane. The least squares fit to data of a linear 
function r = ah gives a = 0.75 ± 0.03 kPa/cm after the col- 
lapse. The shear strength of the unperturbed suspension/soil 
(before the collapse) follows an approximately constant behav- 
ior r ~ 5 kPa until reaching the bottom of the system, which 
was at about 60 cm below the surface for this case. 

The vane rheometer measures both the threshold stress 
and the steady state stress. Here we measure the maximal 
torque necessary to rotate the vane and, knowing its geom- 
etry, it is possible to calculate the threshold stress below 
which the vane does not move [20] . As shown in Fig. [5] be- 
fore destroying the structure, r is essentially constant up 
to the bottom of the basin and then it rapidly increases. 
In order to reproduce this behavior with our model, simu- 
lations with much larger system sizes would be necessary. 
After the system collapsed and the water came out, r lin- 
early increases with depth h: 

r(h) = ah , (1) 

with a = 0.75 ± 0.03 kPa/cm. We conclude from our mea- 
surements that the soil essentially consists of a suspension 
with depth independent shear strength. After the collapse, 
it becomes a soil dominated by the Mohr-Coulomb fric- 
tion criterion for its shear strength. Such a behavior is 
again consistent with our model simulations of a collaps- 
ing soil/suspension. Summarizing, we found a cross-over 
from a yield stress material with a threshold value inde- 
pendent on depth to a Coulomb material after the col- 
lapse. At least, in the experiments the depth-independent 
threshold value can be clearly seen (fig. E|) . 

The analysis of this collapsing soil shows grains with 
a typical size of d — 344 /xm and standard deviation 
<Td — 120 //m as well as the presence of a huge amount of 
cyanobacteria and diatomacea of various types [IB]- Dur- 
ing the drying of the lakes, these organisms are responsi- 



ble for the formation of the elastic and impermeable crust 
which prevents further water from evaporating out of the 
metastable suspension of grains. As reported in previous 
studies [17] , we confirmed that the cohesion force between 
grains inside the suspension is mediated by the cyanobac- 
teria present in the system. 



4 Conclusion 

Our model is capable to simulate the most important fea- 
tures of a real collapsing suspension/soil. Most features 
are important for many types of collapsing suspension/soil 
including wet and dry quicksand, and the model may also 
be applied for the collapse of cohesive powders or snow. 
Shed by bacteria in a highly unstable granular skeleton, 
this suspension/soil can catastrophically collapse. During 
this rapid segregation, objects lighter than water, once 
pushed deep down through the collapsing suspension, can 
be irreversibly trapped under the debris of disassembled 
grains and clusters of grains. According to our simula- 
tion results this behavior appears to be a general feature 
of a collapsing suspension/soil. In our case, the fraction 
of still active bonds with respect to the total number of 
contacts, is a very useful number for describing the col- 
lapsed structures microscopically. Moreover, our model 
results show similar behavior for the shear strength of the 
real material in both unperturbed and perturbed phases, 
showing a cross-over between a yield stress material to a 
Coulomb material after the collapse. As a challenge for 
future work, we intend to investigate the effect of size and 
shape of the intruder on the "drag" force exerted by the 
granular network [21) in the presence of cohesive forces. 
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